function A3 = generate_poisson_matrix(nx0, ny0,alx,aly)
    % 计算x方向的网格间距
    dx = alx / (nx0 + 1);
    
    % 考虑边界，增加网格点数
    ny = ny0 + 2;
    nx = nx0 + 2;
    
    % 定义对称的三对角矩阵 Ax
    Ax = 2*eye(nx0) - diag(ones(nx0-1,1),1) - diag(ones(nx0-1,1),-1);
    Ax = Ax / dx^2; % 转换为差分矩阵算子并除以dx^2
    Ax = sparse(Ax); % 转换成稀疏矩阵

    
    %ny02=ny0/2;
    pi2=2.*pi;

    pi2y=pi2/aly;   %2pi/Ly
    %vky, vkysq: 分别代表y方向的频率值及其平方，系数因子会添加到系数矩阵中
    vkysq=zeros(ny,1);
    vky=zeros(ny,1);
    for j=1:ny0/2+1
        vky(j)=pi2y*(j-1);  %2pi/Ly*(j-1)
        vkysq(j)=vky(j)^2;%(2pi/Ly*(j-1))^2
    end
    for j=ny0/2+2:ny0
        jj=j-ny0;%负指标
        vky(j)=pi2y*(jj-1);
        vkysq(j)=vky(j)^2;
    end

    vkysqz = vkysq;

    % 构建与频率相关的系数向量
    ak2=zeros(ny,1);
    ak2(1:ny0,1)=vkysqz(1:ny0);

    % 计算x方向网格点坐标
    x2=zeros(nx,1);
    for i=1:nx
        x2(i)=dx*(i-1)-alx/2;
    end
    % 结合x和y方向构建大矩阵A4
    A4=kron(ak2(1:ny0,1),eye(nx0))+kron(ones(ny0,1),Ax);
    A4=sparse(A4);
    % 初始化最终矩阵A3
    A3=sparse(nx0*ny0,nx0*ny0);
    % 循环累加构建A3，模拟不同频率模式的叠加
    for i=1:ny0
        % 构建辅助矩阵E和u
        v1=sparse(i,i,1,ny0,ny0);
        E=kron(v1,eye(nx0));
        u=kron([zeros(1,i-1),1,zeros(1,ny0-i)],eye(nx0));
        u=sparse(u);
         % 计算增量矩阵并累加到A3中
        incr=E*A4*u;
        incr=sparse(incr);
        A3=A3+incr;
    end
    % save('A3.mat','A3');
    clear A4;